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Abstract 



A new method for calculating the total energy of Si systems is presented. 
The method is based on the effective-medium theory concept of a reference 
system. Instead of calculating the energy of an atom in the system of in- 
terest a reference system is introduced where the local surroundings are sim- 
ilar. The energy of the reference system can be calculated selfconsistently 
once and for all while the energy difference to the reference system can be 
obtained approximately. We propose to calculate it using the tight-binding 
LMTO scheme with the Atomic-Sphere Approximation(ASA) for the poten- 
tial, and by using the ASA with charge-conserving spheres we are able to treat 
open system without introducing empty spheres. All steps in the calculational 
method is ah initio in the sense that all quantities entering are calculated from 
first principles without any fitting to experiment. A complete and detailed 
description of the method is given together with test calculations of the en- 
ergies of phonons, elastic constants, different structures, surfaces and surface 
reconstructions. We compare the results to calculations using an empirical 
tight-binding scheme. 
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I. INTRODUCTION 



Two parallel developments have changed our possibilities of understanding and predicting 
the energetics and dynamics of condensed-matter systems over the past few years. One is 
the development of ab initio total-energy methods [|l| based on density-functional theory 
0, and the local-density approximation 0. It is now possible to calculate total energies 
and equilibrium configurations of systems with up to a few hundreds of atoms @], and in 
other cases it has been possible to study the molecular dynamics directly During this 

time, there has been a parallel development of approximate and semi-empirical total-energy 
methods |[7|-p!3| which has made it possible to treat the dynamics and thermal properties of 
systems with many thousand atoms with a reasonable accuracy [|T4| , p!5| . With such methods 
it is now possible to treat problems in materials physics where extended defects or long 
range disorder are crucial for the properties that are under study. 

The hope is to develop methods with the accuracy and robustness of the ab initio methods 
which can handle the larger systems on a reasonable time scale. One problem with the ab 
initio methods is the fact that the computer time scales as the cube of the number of atoms. 
New developments of methods where the computer time scales linearly with the size of 



the system have all focussed on localized basis sets [^,0. One important problem here 



is to construct the Hamiltonian. A natural choice would be to use the LMTO method in 
the tight-binding formulation but at the present time these methods usually rely on 



semi-empirical tight-binding models constructed to fit a set of properties. 

In the present paper we introduce a very efficient and accurate method for calculat- 
ing total energies for Si based on density-functional theory. The method is approximate, 
but computationally very efficient (2 orders of magnitude faster than conventional ab initio 
methods). A set of well defined approximations are made in the total energy expression 
where we rely on the variational nature of the generalized total energy functional to com- 
pute energies reliably. Since the approximations are very systematic, we are able to test 
the validity of the assumptions at each stage, and therefore in a controlled way develop a 
hierarchy of models with various levels of approximations. All input terms are calculated 
theoretically, i.e. no fitting to experiment. We have earlier demonstrated the versatility of 
the first levels of approximations in its application to metallic Al [|1^ . 

The method utilizes the effective-medium theory concept of an effective medium or ref- 
erence system. Instead of calculating the total energy of a system of interacting Si atoms 
directly, we associate with each atom in the system a reference system where the energy 
of the atom can be calculated easily and where the surroundings of the atom in question 
are sufficiently similar to those of the system of interest that the energy difference can be 
calculated approximately. Our main approximations are to use a transferable input charge 
density and a transferable effective potential From these approximations the energy 
difference can be calculated with a density- dependent pair potential and an LMTO tight- 
binding Hamiltonian, with all the in-going parameters determined from the transferable 
charge density and effective potential. 

We have made extensive comparisons of the results of this new method with selfconsistent 
calculations of elastic properties, phonons, surfaces, different crystallographic phases, surface 
reconstructions and adatoms on surfaces. In all cases the quality of the results is good, and 
we show that in the cases where the empirical tight-binding method is known to fail, the 
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new method works well. 

The main result of the present paper is the scheme to calculate total energies for Si 
systems, which is summarized in the first part of section IV. The basis of the method is 
the Harris functional and the effective-medium theory concept of a reference system. These 
aspects of the paper are discussed in section II and the first part of section III. In the second 
part of section III we construct the LMTO tight-binding model from which the one-electron 
energy is calculated. The tests of the new method are given in section IV, which also includes 
a discussion of the relation of the present method to the empirical tight-binding method. 



II. THE FIRST LEVEL OF APPROXIMATION: THE OPTIMIZED DENSITY 

AND HARRIS FUNCTIONAL 



A. General remarks 



The so-called Harris functional |^ is a good starting point for investigating the theo- 
retical foundations of the tight-binding method |^3|,24| due to the non-selfconsistent nature 
of the functional and the fact that it only depends on the input charge density. If the input 
density is good enough [^], then this will give a considerable savings in computer time since 
only a single iteration of the Kohn-Sham equations is needed. In this work, we invoke these 
properties of the Harris functional to develop a model potential for Si. 



B. Constructing atomic-like optimized densities 



We have discussed previously [p0| , p6[ the systematic decomposition of a selfconsistent 
total charge density into atomic-like optimized densities 

n(r) = ^^op(r - R^) . (1) 



Here, we use norm- conserving nonlocal pseudopotentials p7| , |28[] within a plane-wave basis 



H at an energy cutoff of 12 Ry to compute the selfconsistent charge densities of bulk Si 
and the ideal (111) surface from which we extract the reciprocal-space components of the 
optimized density Anop{k) (see Appendix I). Reverting to real-space, the optimized density 
shown in Fig. 1 has the well known features of a contraction in the outer core region and 
a sharper attenuation of the tail with some barely discernible Friedel oscillations [BO]. As 



noted before |]2D|, embedding an atom in a homogeneous electron gas at typical metallic 
densities produces very similar features, and this is of interest here because renormalizing 
the atom in this manner is the basis of the effective-medium theory of Jacobsen, Puska and 
N0rskov ig]. 



We have extracted reciprocal-space components of the optimized density for the other 
two principal surface orientations, viz. (100) and (110), and we find that to a good approxi- 
mation, the components fall on the same universal curve. We will now make the assumption 
that Anop{r) is indeed universal and transferable, and in the next subsection we will test this 
ansatz by computing the Harris functional for various different test situations and comparing 
with the selfconsistent results. 



3 



C. Results 



In order to test the accuracy of the approximations at each step in the potential con- 
struction, we will use a data base of test systems. This data base covers silicon in different 
crystallographic structures at the equilibrium volume, phonons, elastic constants, and sur- 
faces. For the surfaces, we have used a supercell geometry with 12 atoms, and for all 
calculations we ensure an adequate sampling |31] of the Brillouin zone. 

In Table 1, we compare the Harris with the selfconsistent results. We also include the 
effective-medium tight-binding (EMTB) results and those of the empirical tight-binding 
model of Ref. which we discuss in the next section. We note that using the Harris 



functional with our choice of input-charge density, constructed from spherically-symmetric 
atomic-like densities (i.e. without bond-charges), there is excellent corroboration with the 
selfconsistent results. Other studies [0,0 using the free atom density as an input into the 
Harris functional have reproduced the bulk selfconsistent results with a similar degree of 
accuracy. 

Our first new result is that the total energies of the ideal surfaces of this semiconducting 
material are accurately derived from a non-selfconsistent calculation. It is interesting to 
note that this is so despite the fact that the exact position of the surface states are known 
to be sensitive to the degree of selfconsistency. The reason is simply that the total energy 
is stationary in the density while the position of the surface states are not. We therefore 
emphasize that the derived potential should only be used for predicting physical quantities 
which are stationary in the density. 



III. THE SECOND LEVEL OF APPROXIMATION: THE EFFECTIVE-MEDIUM 

TIGHT-BINDING MODEL 

A. general remarks 

We will use the effective-medium construction as a basis for making further approxima- 
tions. The effective-medium idea is to first calculate the total energy for a series of reference 
systems, and then for a given system relate each atom to a reference system, and only calcu- 
late the energy difference between the two systems. If the reference system is chosen wisely, 
the energy difference will be small, and can therefore be calculated approximately. To find 
the appropriate reference system, we will introduce the concept of a neutral sphere, defined 
to be a sphere around an atom for which the electron density exactly compensates the pos- 
itive atomic charge. As reference system we will use silicon in the diamond structure with 
a lattice constant such that the neutral-sphere radius is the same in the reference system as 
for the atom in the original system. 

In calculating this energy difference we will utilize that the Hohenberg-Kohn density 
functional can be generalized to a functional E[n,v] which depends on both the density n 



and the potential v |12,23] and which is stationary with respect to independent variations 



of each variable. The general functional can be written as 

E[n, v] = Y, eaH - / n{r)v{r)dr + Eei[n] + E^,[n], (2) 
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where ejyv] denotes the eigenvalues generated by solving the Kohn-Sham equation with the 
potential f , and where -Eeil'^] and -E^cN is the electrostatic and exchange-correlation energy, 
respectively. If the potential is restricted to be a functional of the density the Hohenberg- 
Kohn functional or the Harris functional appear as special cases . 

We have already used the stationary property of the density by calculating the total 
energy using a superposition of atom-based optimized densities. We will show that from this 
approximation the electrostatic and exchange-correlation energy can be transformed into a 
density-dependent pair-potential sum by linearizing the exchange-correlation functional. 

In order to get a simple scheme for calculating the kinetic energy we will use the Atomic- 
Sphere Approximation (ASA). We thereby exploit the stationary property of the the kinetic 
energy functional with respect to variations in the potential by substituting the full potential 
with the spherically-averaged potential within each neutral sphere. Furthermore, we have 
recently shown ||2T| that the spherically-averaged potential of the reference system is very 
similar to that of the real system, and since the kinetic energy functional is stationary in 
the potential it is a good approximation to substitute the spherically-averaged potential 
within each sphere with the potential of the reference system. This last approximation 
transforms the potential-energy contribution to the kinetic-energy difference into a sum of 
density-dependent pair potentials. 

The remaining term is the one-electron energy which we calculate using an LMTO tight- 
binding Hamiltonian. Since the only potential appearing in the problem is now that of 
the reference system, we can precalculate the potential parameters and, by constructing 
an interpolation formula for the structure constants, we obtain a simple density-dependent 
two-center tight-binding Hamiltonian. 

We have now given a brief description of the main approximations used in the construc- 
tion of the Effective-Medium Tight-Binding (EMTB) model. In the following, we give a 
more detailed account of the construction. 



B. The diamond reference system and the neutral-sphere radius 

In the original effective- medium theory of Ref. 0], each atom is viewed as embedded in 
the electron density from the neighboring atoms and, when averaged, this density provides an 
effective medium in the form of a homogeneous electron gas. The role of the homogeneous 
electron gas is to provide a reference system in which the atoms have similar chemical 
surroundings as in the original system. In order to obtain the total energy, corrections due to 
the non smoothness of the charge density have to be included; however these corrections are 
small and can be calculated approximately. We will use the effective-medium construction 
and calculate the total energy E as 

E = Y.e'-^s.) + [E-Y.e^^f{s,)l (3) 

i i 

where we use the neutral-sphere radius Sj of each atom to define the reference system and 
e^^^{s) is the energy of the reference system with neutral-sphere radius s. Note that since 
the electron density is constructed using Eq. (|lD, the choice of identical neutral spheres in 
the two systems is equivalent to the choice of identical embedding density used in Ref. [0]. 
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For silicon in the pseudo-potential scheme the neutral sphere contains 4 electrons, and 
from Eq. (|l|) we obtain the equation 

4 = ^r(d,,,s,), (4) 

3 

where F is the electron-density contribution from atom j to a sphere at site i of radius Sj, 

T{d,s) = jdvAn,^{\Y-d\). (5) 

We solve Eq. @ iteratively for the neutral-sphere radius, using a cutoff distance of Tc = 
11.67ao to terminate the sum (we thereby include five neighbor shells in the equilibrium 
diamond lattice). For this lattice, the nearest-neighbor distance is do = 4.40ao and we 
obtain a neutral-sphere radius of sq = 2.72ao. In Fig. |^ we show the neutral-sphere radius 
of the diamond structure as a function of the nearest-neighbor distance. For comparison, 
we also show the Wigner-Seitz radius of the diamond lattice and the 23 smaller Wigner- 
Seitz radius of a diamond lattice embedded in a bcc lattice with empty spheres. We see 
that in the diamond structure the neutral-sphere radius is substantially smaller than the 
Wigner-Seitz radius. This difference is due to the large regions in the diamond lattice which 
contain almost no charge and therefore are not included in the neutral sphere. In a more 
close packed system like the fee system the neutral-sphere radius is more or less equal to the 
Wigner-Seitz (WS) radius. 

As reference system we will use the diamond structure, and the energy correction there- 
fore vanishes for a diamond lattice or an isolated atom, since the latter can be regarded as 
a diamond lattice with an infinite lattice constant. So by construction the EMTB will give 
the correct cohesive energy, bulk modulus, and equilibrium lattice constant for the diamond 
structure. 

In order to calculate the energy difference to the reference system, we divide the total 
energy into two terms; The kinetic energy (T) and the sum of the electrostatic and exchange- 
correlation energy (G) 

E = Y.^^^f{s,) + [T -Y.^-f{s,)] + [G -Y.r^Si)l (6) 

i i i 

E = J2e''^^{s) + AT + AG. 

i 

So far we have just rearranged the terms in the total energy. The first approximation we 
make is to calculate the kinetic energy in the atomic-sphere approximation (ASA). The 
quality of this approximation depends on how well the full potential can be approximated 
by a superposition of spherically-averaged potentials within atomic spheres. Traditionally, 
the ASA is made using space-filling spheres on the ground that with this choice integrals over 
space may be mapped into integrals over the spheres - the integration of a constant function 
will therefore be correct in the ASA. However, since the density enters all integrals, we will 
use neutral spheres, making the spheres charge conserving instead of volume conserving 
- the integration of a constant function times the density is correct. This choice seems 
more physical since we thereby obtain that both the ASA density and the full density 
contains the correct number of electrons, which is not fulfilled when the ASA is used with 
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volume conserving spheres. In Fig. ^ we show the full selfconsistent potential of silicon in a 
(110) diamond plane, together with the neutral-sphere radius (solid circle) and Wigner-Seitz 
radius (dotted circle). It is clear from the figure that the overlap region of the WS spheres 
penetrates far into the spherically symmetric parts of the potential, which would lead to a 
poor approximation of the full potential. The overlap region of the neutral-spheres, on the 
other hand, very accurately sample only the asymmetric part of the potential and since the 
potentials are superimposed the approximation is enhanced in this region. We will use the 
ASA with neutral spheres, and because we thereby obtain an accurate approximation for 
the full potential we may avoid the traditional usage of empty spheres in ASA calculations 
for open structures. 

We have recently shown that with our choice of reference system the ASA potential of 
a general system is almost identical to the potential of the reference system ^T|, so that 
we may at each site substitute the potential within the atomic sphere with the reference 
potential. Furthermore, due to the variational properties of the energy functional, this only 
leads to errors in the total energy of second order in the potential difference {v — v^'^^). We 
then have for the kinetic-energy difference 

AT^ e.r^VEeieK^O-E / v^'^{r,s,){n{r)-n^^f{r))dr, (7) 

aGocc i i 

= AE,,i + AV, 

where v^^^ is the spherically-averaged potential of the reference system, and eiei{s) is the 
one-electron energy of the reference system with neutral-sphere radius s. 

Let us now summarize the total binding-energy expression in the form used in the 
Effective-Medium theory of Ref . [|I^ . The total binding energy is given by 

E = E, + AEas + AEieu (8) 
= E + [AG + AV] + AEuu 

i 

where i?c, AEas, and AEi^i are called the cohesive function, atomic-sphere correction, and 
one-electron correction, respectively. 

The first term, the cohesive function, is the energy of the reference system, which we 
parametrize using the interpolation formula p4| 



E,(s) =Eo(l+a;)e-", a; = A(s-so). (9) 

In this equation, Eq = —5.83, is the cohesive energy of the equilibrium diamond lattice and 
the parameter, A = 2.047, is determined by the bulk-modulus of the diamond lattice. 

The second term, the atomic-sphere correction, can be viewed as a correction to an 
ASA calculation of the electrostatic and exchange-correlation energy. To see this we use 
the definition of the effective potential v'^^-l' = , where g^'^^^ is the ASA electrostatic and 
exchange-correlation energy of the reference system. We can now identify the — AV^ term of 
Eq. (|^) as the first term in a Taylor expansion of the difference in the ASA electrostatic and 
exchange- correlation energy between the system and the reference systems (AG), and we 
therefore have AEas = AG — AG + 0{[n — n'^'^f]'^) . We see that the atomic-sphere correction 
is the first order correction to a calculation where not only the kinetic energy, but also the 
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exchange and correlation energy have been calculated within the ASA, i.e. the ASA has 
been used for both the potential and the density. In the next section we will show that the 
atomic-sphere correction can be calculated by a density-dependent pair potential. 

The third term, the one-electron correction, is the energy correction due to the difference 
in band structure between the system and the reference system, and in section |1111J| we will 
calculate this term using an LMTO tight-binding model. 

In Table 2 we show the EMTB terms for the test systems of Table 1. The cohesive 
energy and atomic-sphere correction were extracted from a Harris functional calculation 
and the one-electron correction then obtained by subtracting these terms from the total 
energy. With this data base we may not only check the potential by calculating the total 
energy, but we may test the accuracy of each term in the energy separately. For now we will 
just note that the contribution from the cohesive function to the energy of the equilibrium 
structures is very small, indicating that it is the minimum of this function that determines 
the equilibrium volumes. Since the cohesive function depends only on the neutral-sphere 
radius, this implies that all the equilibrium structures have almost the same neutral-sphere 
radius, even though the equilibrium volumes are very different. 



C. Calculating the atomic-sphere correction with a density-dependent pair potential 

We will now decompose the atomic-sphere correction into density-dependent pairwise 
interactions. The interactions depend on the local density through the neutral-sphere radius 
s 

Eas = Y.Vidi^,Si), (10) 
where the pair potential is composed of three parts 

V{d, S) = S) + Vel{d) + V;.,(rf, S), (11) 

which originates from the AV term, the electrostatic energy, and the exchange-correlation 
energy, respectively. From Eq. (|I]) and Eq. (|^ we see that Eq. (^) is exact for the AV^ term 
and the electrostatic energy, with the pair potentials given by 



VMs) = y^r^^(r,s) Ar2,p(|r-d|)(ir, (12) 

V.i{d) = J [vi{r) + i I ^^^dr'] An,,{\r - d\)dr. (13) 

We only include the local part of the pseudo potential vi, since the nonlocal part is canceled 
between the electrostatic energy and the potential part of the kinetic energy. 

To calculate the exchange-correlation energy we use the local-density functional of 
Ref. [35|, and from the decomposition of the density we have 

Exc = Y.j A^op(|r - Ril) S^c{n{v)) dr. (14) 

i 

In order to approximate this term with a pair-potential sum we will have to divide the 
exchange- correlation functional into contributions from each atom, which may be obtained 
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by using a linear expansion for the exchange-correlation functional. We have chosen to 
linearize the exchange-correlation functional around the spherically-averaged density of the 
reference system, n(r, s), and thereby obtain the pair potential 

KcK s) = j An,p{r)-^{n{r, s))An,.^{\r - d\)dr. (15) 

In Table fTT| we compare the exchange-correlation part of the atomic-sphere correction cal- 
culated using this pair potential to an exact evaluation of the exchange-correlation integral. 
We see that the correction generally is underestimated by 20 percent. This suggests that a 
better approximation might be obtained if the exchange- correlation functional is linearized 
at a lower density than the mean density, possibly because it varies more rapidly in the low 
density regime. We have not addressed this problem further, being satisfied with the fact 
that when a factor is allowed for rescaling the pair potential we obtain a good description 
of the exchange-correlation energy, as shown in the third column of Table fTT[ 

We conclude this section by showing in Fig. | the distance dependence of the pair poten- 
tials, and in Fig. ^ the dependence upon the neutral-sphere radius. We see that the distance 
dependence of the sum is nearly exponential even though the distance dependence of the 
individual components is not. The dependence upon the neutral-sphere radius is dominated 
by the contribution from the exchange-correlation part, which approximately scales as s^, 
originating from the scaling of the mean density in the reference system. 



D. Calculating the one-electron correction with an LMTO tight-binding model 

Returning to the expression for the total energy, Eq. (^), we have already given expres- 
sions for the first two terms and need only to calculate the one-electron correction to obtain 
the total energy. This is the most time consuming step in a total-energy calculation since it 
involves the diagonalization of a Hamiltonian. The key number in this context is the num- 
ber of basis functions (N) because the computer time used in conventional diagonalization 
schemes scales as 0(A^'^)[|. A plane-wave basis set is in this respect not very efficient, since 
many plane-waves are needed in order to get a good description of the regions around the 
atomic positions where the electron density varies rapidly. Instead, we will use a partial 
wave method in which the basis functions are augmented with the the local solution of the 
Schrodinger equation around each atom, and therefore only a small basis set is needed. 

However, the partial waves depend on the energy at which the Schrodinger equation is 
solved, and the Hamiltonian thereby becomes energy dependent. This problem is solved very 
elegantly in the linearized band-structure methods where an energy independent basis set is 
obtained by linearizing the solutions of the Schrodinger equation around a fixed energy e^,. 



We will use the Linearized Muffin Tin Orbital (LMTO) method ||T8|, and with these basis 



functions, called LMTO's, the eigenvalues become correct to first order in the difference with 



^In the Car-Parrinello method the scaling is O(A^M^) where M is the system size. However, 
the prefactor is large such that this method generally is two orders of magnitude slower than 
tight-binding methods. 
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the energy e^^ and a systematic expansion exists for the higher order corrections. For now 
we will only consider the first order approximation, and we may then neglect the overlap of 
the orbitals, since the overlap enters as a higher order correction. 

Since we use the ASA for calculating the kinetic energy the first order LMTO Hamiltonian 
becomes especially simple. It separates into potential parameters A", determined from 
the potentials in the atomic spheres, and structure constants 5" which only depend on the 
positions of the atomic spheres 

HTl.v = CtAL.v + (A,^)^5r^^^.^,(AJ,)^ , (16) 

where we use the notation L = Im for the angular-momentum quantum numbers, and we 
will use an sp^ basis set. 

The index a in Eq. (^) denotes the representation we use for the LMTO's. The conven- 
tional LMTO's are the a = representation in which the structure constants have simple 
two center forms 

'S'iL',iL = 5 'S'{gs(T,spcr,ppo-,pp7r} = ^ ^2\/?>x , 12x ^, — 6a; ^} , (17) 

where a; is a relative distance measure given by 

, , d[(si + Si)/2] 
X = dij/Wij , Wij = So -r^^ . (18j 

"0 

In this equation d[s] is a function that returns the nearest-neighbor distance in the diamond 
lattice with neutral-sphere radius s, such that we have the same relative distance Xq = c/q/sq, 
for all diamond lattices. Note that the LMTO Hamiltonian will not depend on the choice 
of w as defined in Eq. (|18D , since the w dependence of the structure constants is canceled 
by a similar term in the potential parameters. 

It is possible to shift to a new representation where the neighboring sites, through a 
screening "charge" a, are used to localize the structure constants |1^. The structure con- 
stants now depend on the local structure through a matrix equation (the LMTO Dysons 
equation) 

^iLJL' = ^iLJL' + X! ^iL,kL"(^l"^kL",jL'- (19) 
kL" 

We use the screening constants a^^p = {0.3072, 0.0316} which are related to the sp-screening 



parameters of Ref. [0 through a scaling of (1.07) . 

With this choice of screening constants we have calculated the structure constants for 
all the test systems of Table 1 using the LMTO Dysons equation, which involves inverting 
a 200x200 matrix. In Fig. |^ we show the Slater-Koster components of each structure 
matrix, as a function of the relative distance measure(x) defined in Eq. (|18|). It is surprising 
how well the structure constants for such different surroundings as surfaces, phonons, and 
different crystal structures, all fall on the same curves. In Ref. interpolation formulas 
for the structure constants of close packed structures were found by using a relative distance 
measure obtained by scaling the distances with the Wigner-Seitz radius. The neutral-sphere 
provides a natural measure, which makes it possible to extend this idea to more open struc- 
tures and surfaces. 
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With our choice of relative distance measure we may now use the data base to construct 
interpolation formulas and thereby obtain a simple and fast scheme for calculating the 
structure constants of a general system. For the ssa and ppw elements we only need to 
include nearest-neighbor elements, while the spa and ppa elements are longer range, and we 
therefore have to include second- nearest neighbor elements. For the interpolation we have 
used the functional form 

S'^ix) = fix) - fix,) -ix- x,)f'ix,), (20) 
/(x)=A(l + A(x-l))e-^^ 

where the first equation ensure that the structure constants go continuous differentiable to 
zero at the cutoff Xc- The relative distance measure, x, is defined in Eq. ( pISD and A, X are 
parameters to be determined from the data base. We fix the parameter A from the nearest- 
neighbor structure constant of the diamond lattice, and determine A by a least-squares fit 
to the nearest-neighbor data points of Fig. |^. The resulting approximations are shown as 
solid lines in Fig. ^ and the values of the parameters are given in Table |rV|. 

The on-site elements of the screened structure matrixes are nonzero. We have from 
Eq. (0) 



CTr iri— > .S^L,kL"^l"^kL",iL>^ (21) 



kL" 



where a" is the on-site element calculated using the approximate structure matrix S*". 
For the diamond structure the on-site element is diagonal and we have a^pidiamond) = 
{1.71, 1.46}, while it will contain off-diagonal components in a general system. However, the 
use of the approximate off-site structure matrixes(S'") in Eq. (|2ll) gives a significant error for 
these components, i.e. the off-diagonal components of cTq are non symmetric and generally 
to large. We will therefore use the following approximation for the on-site structure matrix 

S^,L' = \Kl,l' + ^w,L + 2<,,^,(rfzamonrf)], (22) 

which is a simple average of the symmetrized value of cr" between the system and the 
reference system. 

We now return to the calculation of the potential parameters. These are given by the 
solutions at the linearization energy e^, of the radial Schrodinger equation within each atomic 
sphere [Q. Since we get the potential from the reference system, we only have to calculate 
the potential parameters for the reference system once and for all, and then use the neutral- 
sphere radius to find the potential parameters for a general system. In Fig. |^ and Fig. |^ we 
show the value of the potential parameters as a function of the neutral-sphere radius of the 
reference system. For each system we have chosen at the center of gravity of the occupied 
bands. These data are accurately approximated by the interpolation formulas 



A^(s) = A^(so)e"^-^^°('^~'°) , (23) 
A«(s) = A^(so)e-°-^'^(^-^°) , 
Cp-Cs = 4.67 + (C"(so) - C(so) - 4.67)e-°-^^(^-^«) , 
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where A°(so), A" (sq), (sq) and Cp(so) (unit eV) are the potential parameters in the 
equihbrium diamond lattice listed in Table 0, and s is the neutral-sphere (unit Qq). In 
Fig. 1^ and Fig. ^ the interpolation formulas for the potential parameters are shown as solid 
lines. 

We have now constructed the LMTO tight-binding Hamiltonian directly from the data of 
the ab initio plane-wave calculation. However, due to the small basis set and the incomplete 
description of the interstitial region the calculated bandstructure for the equilibrium dia- 
mond structure does not agree completely with that of the plane-wave calculation shown in 
Fig. ^ (SC). We find the occupied band to be 15 percent to wide and the band gap at the F 
point to be to small. In order to improve the Hamiltonian we have made a least-squares fit of 
the potential parameters to the three lowest eigenstates at the F point and the two lowest at 
the X point. By this procedure we include the effect of the neglected orbitals in an indirect 
fashion. The resulting potential parameters are shown in the second row of Table 0, and 
we see that A"(so) is unchanged while both Ap(so) and (sq) — (sq) have been rescaled 
with a factor 0.77. The corresponding bandstructure is shown in Fig. ^ (EMTB), and we 
now have a good description of the occupied parts of the bands, and the bandgap at the F 
point. 

We are now ready to calculate the one-electron correction (Ai^iei) of Eq. (^. Instead 
of calculating the band energy we will calculate the bond energy [0, because we thereby 
remove any first order dependencies on the onsite elements 

Eiei= E e^-E^^^- (24) 

kdocc i 

In this equation Ni is the site projected occupation and £i the site projected onsite element. 
The bond energy depends only weekly on the shift in the average on-site elements. For 
instance for the three principal surfaces, the largest shift in the average on-site element is 
at the (100) surface, where it is 0.84 eV higher, than in the bulk. Such a shift lowers the 
bond energy 0.19 eV compared to a calculation where the average on-site element is fixed at 
the bulk value. We have chosen to fix the average on-site element to be zero, £i = 0, since 



by this approximation the second term in Eq. |2J vanishes, and this greatly simplifies the 
calculation of forces. Due to this approximation the total energy for the surfaces will be a 
httle to high with the EMTB model. 

In Table |V| we show the one-electron correction obtained with this scheme compared 
to the value of the plane-wave Harris functional calculation. We see that the correction 
generally is overestimated by 2 percent. We can improve the values slightly by scaling all 
potential parameters with this factor, see Table 0, and we then obtain the one-electron 
correction in the last column of Table |V|. We see that for all test systems the accuracy 
is acceptable when compared with the value of the total energy of Table |I|. In the next 
section we will first sum up the ingredients of the EMTB, and then compare it to empirical 
tight-binding schemes by applying both models to various test systems. 



IV. APPLICATIONS 
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A. The EMTB potential 



An EMTB calculation starts with loading precalculated values for the functions 
T{d, s),V{d, s) defined in Eq. (|^) and Eq. (pA]), and the atomic-sphere, e'^J'/{s), and one- 
electron energy term, e[li (s), of the reference system. Next we calculate the neutral-sphere 
radius of each atom, Sj, from Eq. (|). The total energy is given by 

E = Y1 ^c(Si) + Eas-Yl (^as{Si) + E^l " (^lel{Si). (25) 
i i i 

The cohesive function ec(s) is defined in Eq. (H), and the atomic-sphere energy Eas defined 
in Eq. (JT^). The one-electron energy is calculated from the Hamiltonian 

HTl.v = Cri{s.)S,L,L' + (Ar(.,))^^i,L'(^..)(A^(^.))^ , (26) 



where the off-site structure constants are given by the interpolation formula in Eq. ( pOD with 
the parameters of Table and the relative distance measure (x) defined in Eq. ([T8|). The 
on-site elements are given by Eq. ( p2| , pT| , p!7D , using the screening as,p = {0.3072,0.0316}. 



The potential parameters are calculated from the interpolation formulas of Eq. (0) with 
the constraint that the average onsite element at each site is zero {Si = 0), and using the 
values in the last row of Table for A"(so), C"(so). With this model we have calculated 
the total energy for the test systems of Table | and the corresponding results are shown in 
the third column. 



B. Empirical Tight-Binding (EmpTB) 

In a conventional Empirical Tight-Binding scheme the energy function consists of an 
attractive band structure term describing the bonding in the system and a repulsive pair 
potential usually interpreted as arising from the overlap interaction. For the comparison we 



will use the tight-binding model of Goodwin et. al [^, with a fixed cutoff as in Ref. 
and in the following we will denote this model EmpTB. The EmpTB is based on the nearest- 
neighbor tight-binding model of Harrison ^1| , in which the strength of the hopping integrals 



is obtained by fitting the bandstructure, and the level splitting is taken to be identical to that 
of the atom. Harrison assumed an universal decay of the pair potential and hopping integrals, 
such that the only parameter remaining to be determined is the strength of the pair potential, 
which was fixed to give the correct equilibrium lattice constant. The resulting model gives 
an excellent description of the elastic properties of silicon in the diamond structure, however 
the model fails to predict the energies of different silicon phases. Goodwin et. al made 
the model transferable to other structures by introducing a scaling of the hopping integrals, 
and adjusting the level splitting. Their model has four adjustable parameters which were 
fixed to give the cohesive energy and bulk modulus of the diamond and fee structure. This 
scheme has proven very successful for describing systems far from the structures where 



the tight-binding parameters were fitted, and currently most empirical tight-binding schemes 
use a similar functional form. 

The last column in Table | shows the total energies of the test systems obtained with the 
EmpTB. We see that the elastic properties are described very well by the EmpTB, while 
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the total energies of the structures and surfaces are not too accurate. However, note that 
the energy of the crystal structures were calculated at the equilibrium volume as obtained 
from the selfconsistent calculation (in Fig. |10| we show the full phase diagram). 



C. Comparison between the EMTB and EmpTB model. 

In Table | we show for the two models the hopping integrals and sp level splitting in 
the equilibrium diamond lattice. For the EMTB we show the scaled parameters, and for 
the EmpTB we have included the value of ep — used in Ref. [^. In Fig. |^ we show 



the bandstructure of the two models compared to a selfconsistent calculation. For the 
EMTB model the description of the occupied bands is excellent. In the EmpTB model the 
description is reasonable, with largest error for the lowest state at the F point. Also the 
bandgap at the F point is too small in the EmpTB model, this is of importance when the 
model is used together with the 0(N) method of Ref. |T^JT^, where a large bandgap is 
needed in order to make the scheme efficient. 

The difference between the two tight-binding models become apparent when we look at 
the distance dependence of the matrix elements. In the EMTB the distance dependence is 
divided into two parts, scaling of the potential parameters and the structure matrix. When 
we have a uniform compression, without a structural change, only the potential parameters 
change. These scale approximately as [sq/s)^, which is similar to the scaling of Ref. |^ . 
When there is a structural change the scaling should be stronger, and this scaling enters 
through the structure constants. In the EmpTB these two basic scalings are mixed into one 
function. 



In the last column of Table we show the value of the one-electron correction as 
obtained when the EmpTB is used to describe the one-electron energy for the system and 
reference system. Clearly the EmpTB does not describe the one-electron correction, with 
largest discrepancies for surfaces and structural energies. For these systems we have found 
crystal field terms to be important, i.e. off diagonal on-site elements and shifts of the level 
spacings. Such effects are not included in the EmpTB, but enter in the EMTB through 
Eq. (^. 

Besides the one-electron term there is a pair-potential term in both models. However the 
EMTB pair potential is negative, while the EmpTB pair potential is positive. The EmpTB 
pair potential can therefore not only describe the electrostatic and exchange-correlation 



energy, but must include some terms from the kinetic energy. In the work of Harrison |?I 
the pair potential is viewed as an overlap interaction, i.e. it is mainly due to the one-electron 
energy. 

In the EMTB there is an additional term, the cohesive function. The fact that this term 
is not equal to the sum of the reference pair-potential and one-electron energies, illustrates 
the idea behind the reference system: Because the reference system is chosen for a given 
atom so that the environment of the atom is similar to the environment in the real system 
it is possible to calculate the energy difference with a rather crude tight-binding model. 
The largest error is in the shift of the average on-site term of the Hamiltonian, but since 
the potentials are identical in the system and reference system this error for the system 
is canceled by a similar error in the one-electron energy of the reference system, and the 
correct binding energy is then obtained through the cohesive function. 



14 



In the next section we will calculate the phase diagram for silicon and in section [IV E| , [IV F 
we will investigate some of the reconstructions of the (100) and (111) surfaces. 



D. Phase diagram 



In Fig. |T0| we show the total energies versus volume for the diamond, clathrate II 
/3-tin, simple cubic, bcc and fee phases of silicon calculated with the EMTB and EmpTB 
schemes. The EMTB potential gives a good description of the cohesive energy and the 
equilibrium volume for all the phases investigated, while the EmpTB predicts the correct 
energy ordering of all the phases, excluding the clathrate II structure which is lower than 
the diamond structure, but the equilibrium volumes are shifted. The small cutoff of the 
interactions in the EmpTB model becomes visible for the /3-tin and fee phase, i.e. second 
nearest neighbors enter within the cutoff range. In the EMTB scheme we also use a fixed 
cutoff, however in this case we use a relative distance measure which is scale-invariant, and 
this ensures that for a given structure we include the same number of neighbors independent 
of the lattice constant. For the EMTB model we have found crystal field effects to give an 
important contribution to the structural energies. For instance in the fee phase — = 
5.39eV^ compared to — = 5.87eV in the diamond phase. 



E. The (100) surface 

There has been a large effort to understand the different reconstructions of the (100) 
surface, involving a wide range of experimental and theoretical tools since different recon- 
structions occur at different length scales. However, the reconstructions are so complex that 
there are still a lot of unsolved problems, especially for systems involving too few atoms in 
order for elasticity theory to be correct, and too many atoms to be feasible for plane- wave 
methods. 

There is a general consensus that the main building block for all the reconstructions is 



the buckled 1x2 dimer reconstruction. This structure was first predicted by Chadi from 



a tight-binding calculation, and lately also selfconsistent plane-wave calculations p^ , [45[| have 
found this structure to be lower in energy than a symmetric dimer reconstruction. So for a 
potential to give a detailed description of the various reconstructions on the (100) surface, 
it should at least predict the correct 1x2 dimer reconstruction. 



In Table [Vlll[ we show the result of a relaxation of the (100) 1x1 and 2x1 reconstruc- 



tions using the EMTB and EmpTB compared to the results of selfconsistent plane-wave 
calculations. For both models the description of the structure of the 1x2 reconstruction is 
good, but the EMTB gives a too high relaxation energy, while the the relaxation energy of 
the EmpTB is far too small. For the EMTB the driving force for the reconstruction is the 
atomic-sphere correction, and not the one-electron correction. The atomic-sphere correction 
always drives the system to more close packed structures, and in this case that can be done 
without an increase in the one-electron energy. Both models predict an outward relaxation 
of the 1x1 cell, while the selfconsistent calculation predicts an inward relaxation. 
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F. The (111) surface 



Also the (111) surface has many interesting reconstructions, the most famous being the 
7x7 Tagayanagi reconstruction [|^. This reconstruction is built from adatom geometries 
where the adatoms sit on top sites, and studies using selfconsistent plane-wave calculations 
have confirmed the stability of adatoms in top site positions p7| , ^ . For a potential to 
predict the correct reconstructions of the (111) surface it therefore has to describe these 
adatom geometries properly. 

In Table [VIII| we show the result of a relaxation of the (111) 1x1 surface, and of the 
top and hollow site adatom geometries in a \/3X\/3 cell. The selfconsistent numbers are 
from two different references; The top site geometry were calculated in Ref. with a ten 
layer slab and a 12 Ry cutoff. In the other reference the energy of the adatom in both 
the top and hollow positions where calculated but with an eight layer slab using a 6 Ry 
cutoff. We see that both tight-binding models predict the top site to be more stable than 
the hollow site geometry, however the formation energies obtained with the EmpTB are 
far too small, while the EMTB is in excellent agreement with the data of Ref. The 
formation energy obtained with the empirical tight-binding scheme is too low to stabilize the 
adatom, and therefore additional parameters have to be introduced in order to describe the 
7x7 reconstruction ||9|. In Fig. 11 we show the geometry of the adatom in the top position, 
and in Table |VIII| we compare the relaxed coordinates. 



V. SUMMARY AND CONCLUSIONS 

In the present paper we have presented a new method for total energy calculations for 
Si systems. We have discussed in detail the hierarchy of approximations behind the present 
formulation. Starting from a fully selfconsistent solution of the Kohn-Sham equations the 
first level of approximation is to use the Harris functional with transferable, optimized den- 
sities centered at each atomic position deduced from independent selfconsistent calculations 
for different bulk and surface structures. At this level the errors introduced compared to 
the fully selfconsistent calculations are very small. 

The next level of approximation involves the introduction of a reference system. We 
have shown that choosing a reference system with the same neutral-sphere radius (or average 
electron density) as in the real system gives one-electron potentials and potential parameters 
for the LMTO Hamiltonian that can be transferred from the reference system to the real 
system with very little error. The energy difference between the reference system and the real 
system can be calculated from a difference between a density-dependent pair-potential sum in 
the real and the reference system and a one-electron energy difference. The former describes 
to a very good approximation the difference in the electrostatic, exchange correlation and 
part of the kinetic energy. The one-electron energy difference taking care of the rest of the 
energy difference is evaluated using an LMTO tight-binding Hamiltonian. 

The results of the new method are very encouraging. The computational effort is similar 
to empirical tight-binding methods but the results seem to be better. The time consuming 
part of the calculation is the diagonalization of the LMTO Hamiltonian and for this part the 
newly proposed "order N" methods can be used |l6| , p!7| , since the bandgap at the F point is 
well described. 
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APPENDIX A: PARAMETERIZATION OF THE OPTIMIZED DENSITY 

Auopik) ^ Auatomik) + aik'^{k - a2){k - a^) x exp{aik - a^k'^), (Al) 

where ai — 0.078, 02 = 2.584, 03 = 1.465, 04 = 1.969, — 0.962 and Uatom is the atomic 
density. 



17 



REFERENCES 



[1] R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985). 

[2] P. Hohenberg and W. Kohn, Phys. Rev. 136, 3864 (1964), W. Kohn and L.J. Sham, 

Phys. Rev. 140 (4A), A1133 (1965). 
[3] For a review, see e.g. Theory of the Inhomogeneous Electron Gas, S. Lundqvist and N. 

H. March (Eds.) (Plenum Press, New York, 1983). 
[4] I. Stich, M. C. Payne, R. D. King-Smith, J-S. Lin, and L. J. Clarke, Phys. Rev. Lett. 

68, 1351 (1992). K. D. Brommer, M. Needels, B. E. Larson, and J. D. Joannopoulos, 

Phys. Rev. Lett. 68, 1355 (1992). 
[5] A. De Vita, I. Stich, M. J. Gillan, M. C. Payne, and L. J. Clarke, Phys. Rev. Lett. 71, 

1276 (1993). 

[6] P. Margl, K. Schwartz and P. E. Blochl, J. Chem. Phys. 

[7] J. K. N0rskov and N. D. Lang, Phys. Rev. B21, 2131 (1980); J. K. N0rskov ibid. B26, 
2875 (1982). 

[8] M. J. Stott and E. Zaremba, Phys. Rev. B22, 1564 (1980). 

[9] M.S. Daw and M.I. Baskes, Phys. Rev. Lett. 50, 1285 (1983). 
[10] M.W. Finnis and J.E. Sinclair, Phil. Mag. A 50 45 (1984). 
[11] F. Ercolessi, E. Tosatti, and M. Parrinello, Phys. Rev. Lett. 57 , 719 (1986). 
[12] K.W. Jacobsen, J.K. N0rskov and M.J. Puska, Phys. Rev. B35 (14), 7423 (1987). 
[13] S. B. Sinnott, M. S. Stave, T. J. Raeker, and A. E. DePristo, Phys. Rev. B44, 8927 

(1991) . 

[14] O. H. Nielsen, J. P. Sethna, P. Stoltze, K. W. Jacobsen, and J. K. N0 skov, Europhysics 
Letters, in press. 

[15] W. D. Luedtke and U. Landman, Comp. Mater. Scien. 1, 1 (1992). 
[16] M. S. Daw, Phys. Rev. B47, 10895. (1993). 

[17] P.-X. Li, R. W. Nunes, and D. Vanderbildt, Phys. Rev. B47, 1089 (1993). 

[18] O.K. Andersen and O. Jepsen, Phys. Rev. Lett. 53, 2571 (1984). A more detailed 

description of the LMTO method can be found in O.K. Andersen, O. Jepsen, and M. 

Sob, in Electronic Bandstructure and its applications, edited by M. Yussouff (Springer 

Lecture Notes, 1987). 

[19] N. Chetty, K. Stokbro, K. W. Jacobsen, and J. K. N0rskov, Phys. Rev. B 46, 3798 

(1992) . K. Stokbro, N. Chetty, K. W. Jacobsen, and J. K. N0rskov, Proceedings of the 
15*'* Tanigushi symposium, (Springer 1993). 

[20] N. Chetty, K.W. Jacobsen and J.K. N0rskov, Lett. J. Phys., Condens. Matter 3, 5437 
(1991). 

[21] K. Stokbro, N. Chetty, K. W. Jacobsen and J. K. N0rskov, submitted for j. of Physics 
C. 

[22] J. Harris, Phys. Rev. B31 (4), 1770 (1985). 

[23] W.M.C. Foulkes and R. Haydock, Phys. Rev. B39 (17), 12520 (1989). 
[24] O. F. Sankey and D. J. Niklewski, Phys. Rev. B40 , 3979 (1989). 
[25] M.W. Finnis, J. Phys. Condens. Matter 2, 331 (1990). 

[26] l.J. Robertson, M.C. Payne and V. Heine, J. Phys., Condens. Matter 3, 8351 (1991). 
[27] E.L. Shirley, D.C. Allan, R.M. Martin and J.D. Joannopoulous, Phys. Rev. B40, 3652 

(1989). 

[28] D. Vanderbilt, Phys. Rev. B32 (12), 8412 (1985). 



18 



[29] J. Ihm, A. Zunger and M.L. Cohen, J. Phys. C12, 3792 (1979). 

[30] M.J. Puska, R.M. Nieminen, M. Manninen, Phys. Rev. B24, 3037 (1981). 

[31] H.J. Monkhorst and J.D. Pack, Phys. Rev. B13, 5188 (1976). 

[32] A.J. Read and R.J. Needs, J. Phys. Condens. Matter 1, 7565 (1989). 

[33] H.M. Polatoglou and M. Methfessel, Phys. Rev B37 (17), 10403 (1988). 

[34] J.H. Rose, J.R. Schmidt, F. Guinea, and J. Ferrante, Phys. Rev. B 29, 2963 (1984). 

[35] D.M. Ccperlcy and B.J. Alder, Phys. Rev. Lett. 45, 566 (1980). J. Perdew and A. 

Zunger, Phys Rev B23, 5048 (1981). 
[36] J.C. Slater and G.F. Koster, Phys. Rev. 94, 1498 (1954). 

[37] A. P. Sutton, M. W. Finnis, D. G. Pettifor and Y Ohta. J. Phys. C 21, 35 (1988). 

[38] Mvffin Tin Orbitals and Electronic StrucfMre by H.L. Skrivcr, Springer- Verlag (1984). 
[39] L. Goodwin, A.J. Skinner, and D.G. Pettifor. Europhys. Lett. 9, 701 (1989). 
[40] G.Z. Wang, G.T. Ghan, and K.M. Ho, Phys. Rev. B45, 12227 (1992). 
[41] W. A. Harrison, Phys. Rev. B27, 3592 (1983). 

[42] The clathrate H structure is a super-cell geometry consisting of 34 atoms all sp^ bonded, 
which has been proposed as a new carbon phase by R. Nesper, K. Vogel, P.E. Blochl 
Angew. Ghem. 32, 701 (1993). We were made aware of the failure of the EmpTB 
model for the Glathrate II structure by A. A. Demkov , and O. F. Sankey, Private 
communications, who also provided us with the total energy of the structure from a 
plane- wave self consistent calculation. 

[43] D.J. Ghadi, Phys. Rev. Lett. 43 (1), 43 (1979). 

[44] N. Roberts and R.J. Needs, Surface Science 236, 112 (1990). 

[45] J. Dabrowski and M. Scheffler, App. Surf. Sci. 56, 15 (1992). 

[46] K. Takayanagi, Y. Tanishiro, M. Takahashi, and S. Takahashi, J. Vac. Sci. Technol. A3, 
1502 (1985). 

[47] J. E. Northrup Phys. Rev. Lett. 57, 154 (1986). Note that this calculation was performed 
using a smaller unit cell than ours, and with a 6 Ry plane-wave cutoff and k-point 

sampling. 

[48] R. D. Meade and D. Vanderbildt, Phys. Rev. B40, 3905 (1989). The authors have not 
calculated the total energy of the ideal (111) surface, so in order to obtain energies rela- 
tive to the ideal surface we use the total energy from our own selfconsistent calculation 
of the clean (111) surface. 

[49] D. J. Ghadi, Phys. Rev. B 30, 4470 (1984). 

[50] M.T. Yin and M.L. Gohen, Phys. Rev. B24, 2303 (1981). 



19 



FIGURES 



FIG. 1. The free atom density (solid) and optimized density (dashed) for Si. 

FIG. 2. The neutral-sphere radius (solid) for silicon in the diamond lattice as a function of the 
nearest neighbor distance. The dashed lines show the Wigner-Seitz radius of the diamond lattice, 
and of an inscribed bcc lattice. 

FIG. 3. The selfconsistent potential of silicon in the diamond lattice projected onto the (110) 
plane. The solid circle show the radius of the neutral-sphere radius, while the dashed circle show 
the Wigner-Seitz radius. 

FIG. 4. The density-dependent pair potential, and its three components, used for calculating 
the atomic-sphere correction as function of distance, with the neutral-sphere radius fixed at sq. 
The inlet shows the logarithm of the pair potential. 

FIG. 5. The pair potential as a function of the neutral-sphere radius (s), with the distance 
fixed at (Iq. 

FIG. 6. The crosses show the Slater-Koster components of the structure constants for the 
structures in Table | as a function of the relative distance measure, x, defined in Eq. (^). The 
structure constants were calculated using Eq. ([T9|). The tick marks on the horizontal axis indicate 
the relative distance in the diamond, simple cubic, fee and bcc structure. The solid lines show the 
value of the structure constant as obtained from the interpolation formula of Eq. (|20|) . 



FIG. 7. The dots show the calculated potential parameters A" for the diamond reference 
system, the solid line is the approximation obtained with the interpolation formulas of Eq. (p4|). 



FIG. 8. The dots show the potential parameters C", and the solid line is the approximation 
obtained with the interpolation formulas of Eq. (|24 



FIG. 9. The figure show the bandstructure of silicon calculated with three different methods. 
Starting from the top the calculations are: Self-consistent plane- wave calculation (SC), the EMTB 
model and the Empirical Tight-Binding model(EmpTB) of Ref. [39|. 



FIG. 10. The triangles in both figures show the energies of the diamond, clathrate II, /3-tin, 
simple cubic, bcc and fee structure, in that order, at their respective equilibrium volumes, calculated 
selfconsistently. The energy of the clathrate II structure is from Ref. |^^. The solid lines in the 
upper figure show the energies of the structures calculated using the EMTB model, while the 
energies in the lower figure where calculated using the empirical tight-binding model of Ref. |39|. 



The cusps on the curves in the lower figures are caused by second-nearest neighbors entering within 
the cutoff distance. 
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FIG. 11. The geometry of the top site adatom on the (111) surface. 
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TABLES 

TABLE L The selfconsistent (SC), the Harris functional (H) and the Effective- Medium 
Tight-Binding (EMTB) results for the lattice constant ao, bulk modulus B and cohesive energy 
Ec for Si in the diamond structure. The energy relative to the diamond phase of the /3-tin, simple 
cubic, bcc and fee structure at the equilibrium lattice constant as obtained from the selfconsistent 
calculation. The phonon frequencies of the transverse acoustics phonon at the X point, the trans- 
verse optical phonon at the X point, the longitudinal acoustic and optical phonon at the X point 
LAO{X), and the longitudinal and transverse optical phonon at the T point LTO{T). The three 
cubic elastic constants. The energies of the ideal principal surface orientations. In the last column 
(EmpTB) we show the values obtained with the empirical tight-binding model of Ref. p9| . 



Quantity 


units 


SC 


H 


EMTB 


EmpTB 


a) DIAMOND BULK 












ao 


A 


5.395 


5.380 


5.380 


5.42 


B 


MBar 


0.99 


0.93 


0.93 


1.04 


Ec 


eV 


5.85 


5.83 


5.83 


4.70 


b) STRUCTURES 












/3-tin (4.76 A) 


eV/atom 


0.26 


0.25 


0.27 


0.60 


sc(2.51 A) 


eV/atom 








u. / 


bcc(3.03 A) 


eV/atom 


0.55 


0.52 


0.73 


1.31 


fcc(3.79 A) 


eV/atom 


0.56 


0.54 


0.89 


1.32 


c) PHONONS 












vta{X) 


THz 


4.3 


4.2 


4.0 


5.0 




THz 


14.1 


13.3 


19.0 


16.6 


vloa{X) 


THz 


12.3 


12.2 


12.7 


14.4 


I^LTOO^) 


THz 


15.7 


15.6 


19.6 


18.5 


d) ELASTIC CONSTANTS 










Cii 


Mbar 


1.70 


1.65 


1.43 


1.46 


Cl2 


Mbar 


0.72 


0.62 


0.90 


0.78 


C44 


Mbar 


1.10 


1.10 


1.69 


1.23 


e) SURFACES 












(100) 


eV /atom 


2.19 


2.10 


2.39 


1.65 


(110) 


eV/atom 


1.27 


1.28 


1.55 


1.43 


(111) 


eV /atom 


1.45 


1.43 


1.56 


1.43 
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TABLE II. The effective-medium components of the Harris functional results of Table |. The 
second column show the lattice constant of the structure and the displacement used for the frozen 
phonon and elastic deformation calculations. 



Quantity 




H 


Ec 




AGxc 


AV 


A^lei 


/3-tin 


(4.76 A) 


0.25 


0.00 


-1.39 


-0.80 


0.67 


1.78 


sc 


(2.51 A) 


0.41 


0.01 


-0.99 


-0.64 


0.37 


1.67 


bcc 


(3.03 A) 


0.52 


0.01 


-1.94 


-1.13 


0.96 


2.62 


fee 


(3.79 A) 


0.54 


0.00 


-1.75 


-1.16 


0.61 


2.84 


TA(X) 


(0.04) 


0.025 


0.002 


-0.053 


-0.022 


0.035 


0.062 


TO(X) 


(0.02) 


0.059 


0.001 


-0.062 


0.023 


0.014 


-0.044 


LOA(X) 


(0.02) 


0.099 


0.029 


-0.050 


0.007 


0.096 


0.017 


LTO(r) 


(-0.01) 


0.145 


0.002 


-0.133 


0.031 


0.277 


-0.032 


LTO(r) 


(0.01) 


0.099 


0.001 


-0.068 


0.025 


0.150 


-0.009 


Cn 


(0.06) 


0.0361 


0.0229 


-0.0060 


-0.0034 


0.0051 


0.0174 


2(Cii — C12) 


(0.06) 


0.0450 


0.0005 


-0.0257 


-0.0117 


0.0170 


0.0650 


<-/44 


(0.06) 


0.0241 


0.0003 


-0.0235 


0.0050 


0.0463 


-0.0039 


(100) 




2.10 


0.40 


0.36 


0.68 


0.32 


0.34 


(110) 




1.28 


0.10 


0.20 


0.52 


0.23 


0.23 


(111) 




1.43 


0.10 


0.24 


0.58 


0.24 


0.27 



TABLE III. The table shows the value of the AGxc term for the structures of Table ||. The 
column denoted H is the value of the term calculated using the Harris functional, the column 
denoted EMTB is the result obtained using the approximations of the EMTB model, and in the 
third column a parameter has been allowed in order to scale all the energies. 



Quantity 


H 


AG,.e 

EMTB 


xl.28 


/3-tin 


(4.76 A) 


-0.80 


-0.62 


-0.80 


sc 


(2.51 A) 


-0.64 


-0.48 


-0.62 


bcc 


(3.03 A) 


-1.13 


-0.90 


-1.16 


fee 


(3.79 A) 


-1.16 


-0.94 


-1.21 


TA(X) 


(0.04) 


-0.022 


-0.015 


-0.020 


TO(X) 


(0.02) 


0.023 


0.016 


0.020 


LOA(X) 


(0.02) 


0.007 


0.001 


0.002 


LTO(r) 


(-0.01) 


0.031 


0.023 


0.029 


LTO(r) 


(0.01) 


0.025 


0.018 


0.023 


Cii 


(0.06) 


-0.0034 


-0.0025 


-0.00339 


2(Cii - C12) 


(0.06) 


-0.0117 


-0.0082 


-0.0105 


C44 


(0.06) 


0.0050 


0.0030 


0.0038 


(100) 




0.68 


0.31 


0.40 


(110) 




0.52 


0.25 


0.32 


(111) 




0.58 


0.26 


0.33 
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TABLE IV. The value of the parameters in Eq. (|20|). The parameter Xc determines the range 
of the Hamiltonian, which for the ssa and ppvr element is nearest neighbor in the diamond lattice, 
while it is second nearest neighbor for the spa and ppa element. The ^"((io/so) is the structure 
constants of the equilibrium diamond lattice as obtained from the LMTO Dysons equation, while 
A is obtained from a least squares fit to the data points of Fig. ^ 



quantity 


SSfJ 


spo" 


ppcr 


ppvr 




2.60 


2.95 


2.95 


2.60 




-0.938 


1.690 


3.279 


-1.025 


A 


2.40 


2.85 


2.76 


4.10 



TABLE V. The LMTO potential parameters for the equilibrium diamond lattice. The values 
in the first row are obtained by solving the radial Schrodinger equation within the atomic-sphere, 
the values in the second row are obtained by a least squares fit to the bandstructure of silicon, and 
the values in the last row are those of the second row but rescaled with 0.98. 



quantity 




A?(5o) 




calc. 


10.73 


1.954 


0.959 


fit 


7.91 


1.954 


0.738 


xO.98 


7.75 


1.915 


0.732 
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TABLE VI. The table shows the value of AEiei terms for the structures of Table The 
column denoted H is the value of the term calculated using the Harris functional, the column 
denoted EMTB is the result obtained using the approximations of the EMTB model, and in the 
third column a parameter has been allowed in order to scale all the energies. The last column 



(EmpTB) shows the values obtained with the empirical tight-binding model of Ref. |39] 



Quantity 


H 


EMTB 


xO.98 


EmpTB 


/3-tin 


(4.76 A) 


1.78 


1.81 


1.77 


0.38 


sc 


(2.51 A) 


1.67 


1.69 


1.66 


0.81 


bcc 


(3.03 A) 


2.62 


2.90 


2.84 


0.29 


fee 


(3.79 A) 


2.84 


3.29 


3.23 


0.44 


TA(X) 


(0.04) 


0.062 


0.058 


0.057 


0.057 


TO(X) 


(0.02) 


-0.044 


0.022 


0.022 


0.049 


LOA(X) 


(0.02) 


0.017 


0.031 


0.031 


0.062 


LTO(r) 


(-0.01) 


-0.032 


0.053 


0.052 


0.1081 


LTO(r) 


(0.01) 


-0.009 


0.055 


0.054 


0.1233 


Cii 


(0.06) 


0.0174 


0.0162 


0.0159 


0.0137 


2(Cii — C12) 


(0.06) 


0.0650 


0.0530 


0.0520 


0.0422 


C44 


(0.06) 


-0.0039 


0.0126 


0.0123 


.0168 


(100) 




0.34 


0.91 


0.89 


1.15 


(110) 




0.23 


0.68 


0.67 


1.23 


(111) 




0.27 


0.65 


0.64 


1.23 



TABLE VII. The tight-binding parameters of the equilibrium silicon structure for the EMTB 
and the Empirical Tight-binding (EmpTB) scheme of Ref. [^. For the EMTB we show the rescaled 
parameters. In Ref. |^9| the four hopping integrals for the EmpTB were taken from Ref. [^l| and 
the level splitting fitted to the fee - diamond energy difference. The number in parenthesis is the 



value of Ref. |4|]. 



model 


Hssa 


Hspa 


Hppcr 


HppiY 




EMTB 


-1.79 


1.99 


2.37 


-0.74 


5.87 


EmpTB 


-1.82 


1.96 


3.06 


-0.87 


8.295(6.83) 
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TABLE VIII. Energies and structures of the (100) and (111) surface obtained from selfconsis- 
tent, EMTB, and EmpTB calculations. The 7 denote the surface energy per 1x1 cell, and A7 the 
surface energy relative to the ideal 1x1 cell. For the relaxed 1x1 geometries we show the relative 
relaxation of the first layer atoms, A12. For the (100) 1x2 structure we show the dimer bond length, 
r^, and buckling angle, 6. For the (111) ^/2>X^/2> T4 structure we show the relaxation toward the 
adatom axis (dotted line in Fig. [Tl| ), 5r, and the relaxation in the vertical direction, 5z. The atom 
numbers refer to Fig. |ll|. 



geometry 


quantity 


unit 


SC 


EMTB 


EmpTB 


THE (100) SURFACE 










Ideal 1x1 


7 


eV 


2.19 


2.39 


1.65 


Rel. 1x1 


A7 


eV 


-0.03" 


-0.02 


-0.03 




A12 


% 


-5.1" 


4.8 


4.5 


Rel. 1x2 


A7 


eV 


-0.85'' 


-1.04 


-0.39 




e 


Degree 


(16°)^ 
4.28^^ 


19° 


I50 




rd 


ao 


4.47 


4.62 


THE (111) 


SURFACE 










Ideal 1x1 


7 


eV 


1.45 


1.56 


1.43 


Rel. 1x1 


A7 


eV 


-0.06^ (-O.l?'^) 


-0.01 


-0.02 




A12 


% 


-27^= 


-7 


3 




A7 


eV 


-0.27^= (-0.28'^) 


-0.27 


-0.04 




rd 


ao 


5.01^= (4.70^) 


4.89 


4.95 




5r{2) 


ao 


-0.35^= (-0.28'^) 


-0.35 


-0.23 




5z{2) 


ao 


0.03^= (-0.15"^) 


0.10 


0.12 




5z{2>a) 


ao 


-0.7V (-0.74^^) 


-0.65 


-0.57 




5z{3b) 


ao 


0.38^ (0.18"^) 


0.18 


0.21 




(5z(4a) 


ao 


-0.54^= (-0.48^^) 


-0.58 


-0.62 




6z{4b) 


ao 


0.23^ (0.11"^) 


0.06 


0.12 




Sr{5) 


ao 


0.10^= 


0.09 


0.13 




5z{5) 


ao 


o.or 


0.03 


-0.04 




A7 


eV 


-0.17'' 


-0.13 


-0.02 



"Reference 501 Reference [45 '^Reference 48 Reference 47 
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0.06 




radial distance (aj 



3.6 - 
3.4 - 



_ ws 

03° 3.2 - 




2 2 ' ' ' ' ' ' ' ' ' ' ' 

■ 4.0 4.2 4.4 4.6 4.8 5.0 

neighbor distance (aj 




(1,1,0) 



0.5 




, I I I I I I I I I I I 

2.5 2.7 2.9 3.1 3.3 3.5 

s (ao) 



-15.0 



-20.0 



> 

3. -25.0 

o 



-30.0 



-35.0 




neutral-sphere radius (aj 



